Certifying the Potential Energy Landscape 
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It is highly desirable for a numerical approximation of a stationary point for a potential energy 
landscape to lie in the quadratic convergence basin of that stationary point. However, it is possible 
that an approximation may lie only in the linear convergence basin, or even in a chaotic region, and 
hence not converge to the actual stationary point when further optimization is attempted. Proving 
that a numerical approximation will quadratically converge to the associated stationary point is 
termed certifying the numerical approximation. We employ Smale's a-theory to stationary points, 
providing a certification that serves as a mathematical proof that the numerical approximation 
does indeed correspond to an actual stationary point, independent of the precision employed. As a 
practical example, employing recently developed certification algorithms, we show how the a-theory 
can be used to certify all the known minima and transition states of Lennard- Jones LJjv atomic 
clusters for N = 7, . . . , 14. 



Introduction: The surface defined by a potential, V(x), 
with x = (xi, . . . ,x n ), is called the potential energy land- 
scape (PEL) of the corresponding physical or chemical sys- 
tem p], 0| . The critical points of a PEL, defined by the 
solutions of the equations dV(x)/dxi = for i = 1, . . . , n, 
provide important information about the landscape. These 
critical points, the stationary points (SPs) of the PEL, can 
be classified according to the number of negative eigenval- 
ues of the Hessian matrix, Hij = d 2 V(x)/dxidxj, evalu- 
ated at the SPs: the SPs with no negative eigenvalues are 
minima, and the SPs with exactly / negative eigenvalues 
are called saddles of index I. SPs at which H has one or 
more zero eigenvalues in addition to those determined by 
translational and rotational symmetry are called singular 
SPs, or non-Morse points. 

Except for rare examples, it is not usually possible to ob- 
tain the SPs analytically, and one has to rely upon comput- 
ing numerical approximations by solving the corresponding 
equations. For a numerical approach, "solve" means "to 
compute a numerical approximation of the associated so- 
lutions." Once a numerical approximation of a solution is 
obtained, it is heuristically validated. Two standard ap- 
proaches are to monitor iterations of Newton's method and 
to substitute the approximations into the equations to see 
if they are satisfied up to a chosen tolerance. Although 
such a validation usually works well in practice, it does 
not guarantee that the numerical approximation will indeed 
converge quadratically to the associated solutions using ar- 
bitrary precision. In other words, even if a numerical ap- 
proximation is heuristically validated, it could correspond 
to a nonsolution when using higher precision. Additionally, 
Newton iterations may have unpredictable behavior, such 
as attracting cycles and chaos, when applied to points that 
are not in a basin of attraction [3-6] of some solution. 

If the given system is a set of polynomial equations, then 
one can use numerical polynomial homotopy continuation 
0-E3] to compute all the isolated solutions. Due to the 
numerical computations used with this method, one ob- 
tains numerical approximations of the isolated solutions 
and hence the aforementioned difficulties also arise. 

A proper validation of a numerical approximation is termed 



certification, i.e., a verification that the given numerical ap- 
proximation will converge quadratically to the nearby as- 
sociated solution using arbitrary precision. Roughly speak- 
ing, quadratic convergence doubles the number of correct 
digits after each iteration, so that the associated solution 
can be approximated to a given accuracy efficiently. Start- 
ing in the 1980's, Smale and others developed a method 
that certifies a numerical approximation as an actual solu- 
tion of the system in the following way [l8|. For a given sys- 
tem of equations / = and a given point x* , one computes 
a number a(f, x*) which, if it less than (13 — 3vl7) /4 « 
0.157671, guarantees that Newton's method starting from 
x* will quadratically converge to a solution of / = 0. More- 
over, by using such a certification scheme, we ensure that 
our numerical approximations of solutions are good enough 
so that more accurate approximations of the solutions can 
be obtained easily and efficiently. 

Smale's a-Theory: We summarize Smale's a-theory fol- 
lowing Ref. [l^, where we restrict ourselves to systems of 
equations that have the same number of equations as vari- 
ables, termed square systems. We should also emphasize 
that Smale's a-theory is usually used to certify complex 
solutions for systems of polynomial equations, so we start 
with the key points of the theory for this case [lH . Since 
more can be said about real solutions, we will discuss them 
separately below, as well as generalizations to other types 
of nonlinear equations, such as the those involving expo- 
nentials and trigonometric functions. 

For a system / of n multivariate polynomial equations in 
n variables, we denote the set of solutions of / = as 
V(/) := {z <E C n \f(z) = 0} and the Jacobian of / at x as 
J/(x). Consider the Newton iteration of / starting at x 
defined by 



x-J/(x)-V(x) 



if J/(x) is invertible, 
otherwise. 



For k > 1, the fc-th Newton iteration is simply 
ATjf(x) :=Nf0...oNf(x). 



k times 
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A point x G C" is called an approximate solution of / with 
associated solution z € V(/) if, for each k > 1, 



Km 



where || • || is the standard Euclidean norm on C". In other 
words, x is an approximate solution to / if it is in the 
quadratic convergence basin defined by Newton's method 
of some solution z. The key to Smale's a-theory, as shown 
in the following theorem, is a sufficient condition for prov- 
ing that a given point is an approximate solution without 
knowledge about z. 

Theorem: If a(f, x) < (13 — 3\/17) /4 for a square poly- 
nomial system / and point x, then x is an approximate 
solution to / where 



W,x) 
7(/,x) 



/3(/,x) 7 (/,x), 
II^Xx)- 1 /^)!!, and 



sup 

k>2 



Jf^)- 1 D k f^) 
15! 



In 7(/, x), the term Z? fe /(x) is the symmetric tensor whose 
components are the partial derivatives of / of order k. 
Additionally, for convenience, if at some x <G V(/) where 
J/(x) is not invertible, then a(f,x) := 0, (3(f, x) := and 
7 (/, x) := oo. If x ^ V(/) such that J/(x) is not invertible, 
then a(f, x), (3(f, x) and 7(/, x) are taken as oo. Finally, if 
x is an approximate solution of /, then ||x— z|| < 2/3(/, x) 
where z e V(/) is the associated solution to x. 

We remark that since this theorem provides a sufficient 
condition for a point to be an approximate solution, the 
set of certifiable approximate solutions is generally much 
smaller than the set of approximate solutions. However, 
for a true approximate solution, a few Newton iterations 
usually generate a point that one is able to certify. 

Given two approximate solutions Xi and X2 , one often needs 
to verify that the corresponding associated solutions zi 
and Z2 are distinct. One way to verify this is by using 
the triangle inequality together with ||xi — Zj|| < 2(3(f, Xj). 

Other Nonlinear Systems: The above theorem was 
actually proved with "polynomial" replaced by "analytic." 
However, we present it in this fashion since, in the poly- 
nomial case, 7 (/, x) is actually defined as a maximum over 
finitely many terms, since only finitely many partial deriva- 
tives can be nonzero. In fact, it can be bounded above 
based on the coefficients of /, the degree of the polyno- 
mials in /, and J/(x). Nonetheless, 7(/, x) can also be 
bounded above for other nonlinear systems, in particular, 
systems of polynomial-exponential equations |20j ] . A sys- 
tem is polynomial-exponential if it is polynomial in both 
the variables x±,...,x n and finitely many exponentials of 
the form e aXi where a € C. Many standard functions such 
as sin(x), cos(x), sinh(x), and cosh(x) can be formulated 
as systems of polynomial-exponential functions since they 
are indeed polynomial functions of e ax for suitable a € C. 

Real Solutions: For a square system / such that Ay- 
defines a real map, i.e., Aj(x) is real whenever x is real, 
then Smale's a-theory can be extended to provide more 
information about real solutions [H|. For potential energy 
landscapes, when the potential energy function l^(x) is real 



for real x, the corresponding Newton iteration is always a 
real map. In this case, one can determine the reality of the 
associated solution y from any approximate solution x. If 
x is real, then y must also be real. However, if x is not 
real, one can show that y is real by showing that x and 
its real part, namely (x + x)/2 where x is the conjugate of 
x, have the same associated solution, namely y. To show 
that y is not real, one simply has to show that x and its 
conjugate x have distinct associated solutions, namely y 
and y, which is shown using ||x — y|| < 2/3(/, x) with the 
triangle inequality. 

We use a recently developed practical implementation of 
the a-theory, called alphaCertified, for certifying solu- 
tions to systems of equations (l9l . [2p| . When using exact 
rational arithmetic, the implementation of a-theory in al- 
phaCertified is rigorous and can be taken as a mathemat- 
ical proof of the computed results. Hence, this approach 
provides an alternative to other analytic or symbolic com- 
putations. The algorithms are also implemented in arbi- 
trary precision floating point arithmetic in alphaCerti- 
fied, which provides certified results up to round-off errors. 

An Illustrative Example: As a demonstration of com- 
puting a(/, x), f3(f,x), and 7 (/, x) for a single coordinate, 
consider the univariate polynomial f{x) = x 4 — 1. In this 
simple case, we can actually compute these quantities as a 
function of a variable x rather than at a specific value. We 
will assume x ^ since f'(x) = 4x 3 is zero if and only if 
x = and /(0) ^ 0. Clearly, /3(f,x) = \x - .t" 3 |/4. Now, 
in the univariate case, the term D k f(x) in j(f, x) is simply 
the fc-th derivative of / at x, i.e., f^(x). Since / has de- 
gree 4, we only need to take the maximum over k = 2, 3, 4 
to compute "/(f,x). One can easily verify that the maxi- 
mum is attained at k = 2 with 7 (/, x) = 3\x~ 1 \/2. Thus, 
a(/, x) = 3|1 — x~ 4 \/8 for any x ^ 0. 

For example, a(/, 2.5) = 0.3654 so that x = 2.5 cannot be 
certified as an approximate solution. In fact, x = 2.5 is 
indeed outside of all of the quadratic convergence basins. 
However, since a(f, 1.1) = 0.11887, x = 1.1 is certifiably 
an approximate solution of / = 0. In this case, we know 
that the associated solution is z = 1 and the following table 
confirms the quadratic convergence for small values of k. 



k 


1 


2 


3 


4 


5 


-log 10 (\\N*(x)-z\\) 


1.89 


3.62 


7.06 


13.94 


27.70 


-log 10 (\\x-z\\/2* h -^ 


1.30 


1.90 


3.11 


5.52 


10.33 



Example With Close Roots: The n-th Chebyshev 
polynomial of the first kind is well-known to have n roots 
between —1 and 1. These roots, called Chebyshev nodes, 
are located at Xi = cos [(2i — l)7r/2n] for i = 1, . . . ,n. We 
can use this example to demonstrate how small pertur- 
bations in a numerical approximation can change which 
root Newton's method will converge to. This chaotic be- 
havior can be avoided by using certification. In partic- 
ular, the following table considers selected values where 
f(x) = cos(50 cos -1 x) is the 50-th Chebyshev polynomial 
of the first kind. 
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X* 


lim N)(x*) 


QQ7 




QQ7Q 


X3 — cos(57r/100) 


QQ7QQ 




0.997999 


xq = cos(ll7r/100) 


0.998001 


x 6 = cos(11tt/100) 


0.99801 


x 9 = cos(17tt/100) 


0.9981 


xi = cos(7r/100) 


0.998 


xq = cos(ll7r/100) 



Miiller-Brown Surface: The Muller-Brown surface [21| 
is a well-known model landscape [22T-r25]|. It is defined as 

4 

V(x, y) = J2 A > ex P ( a i( x ~ x °f +bi(x- x°)(y - y°) + a{y - y°) 2 ) 

i=l 

where 

A = (-200,-100,-170,15), a = (-1,-1,-6.5,0.7), 

b = (0,0,11,0.6), c = (-10,-10,-6.5,0.7), 

x° = (1,0,-0.5,-1), y° = (0,0.5,1.5,1). 

Since W = [dV/dx, dV/dy] involves polynomials as ex- 
ponents, we simply add new variables to produce an equiv- 
alent polynomial-exponential form as 



f (x, J/, 21, ... , Zi, t»l, ... , W4) = 

Etl A z w l (2a z (x - x°) + b,(y - vl)) 
Ei=l A iWi (bi(x - x°) + 2 Ci (y - y°)) 
ai (x - x°) 2 + bi(x - x°)(y - y°) + c t (y - y°) 2 - z u i = 1, . . . ,4 
exp(zi) — Wi , i = 1, . . . , 4 



Given (x,y), we obtain values of z, and u>^ based on the 
last eight functions in / and then try to certify the result. 
In particular, the following table presents five numerical 
approximations of SPs for V along with an upper bound 
on the value of a(f, ■) and an approximation of ■). 







upper bound 


approximation 


X 


y 


of a(f, ■) 


of /?(/,-) 


-0.5582236346 


1.441725842 


0.0140 


4.84 ■ 10" 9 


0.6234994049 


0.02803775853 


0.0460 


1.94 ■ 10~ 9 


0.212486582 


0.2929883251 


0.0437 


3.05 ■ 10~ 9 


-0.8220015587 


0.6243128028 


0.0006 


6.94 ■ lO" 10 


-0.050010823 


0.4666941049 


0.0068 


2.89 ■ 10" 9 



In particular, based on the upper bounds on a(f, •) com- 
puted by alphaCertified, each point is indeed an approx- 
imate solution. The bounds on the distance from each nu- 
merical approximation to the corresponding approximate 
solution based on /3(/, •) show that each one must corre- 
spond to a distinct approximate solution. Thus, we have 
proved that the five numerically approximated SPs are in- 
deed in the quadratic convergence basin of distinct SPs. 
Lennard-Jones Clusters: We now consider one of the 
most studied family of systems in molecular science, namely 
atomic clusters of N atoms bound by the Lennard-Jones po- 
tential [2(|, denoted LJ^r. The pairwise potential between 
interacting particles is defined as 



V N 



N N 

i=i j=i+i 



a 



where e is the pair well depth, 2 1 / 6 cr is the equilibrium pair 
separation, and 

n,j = sj (xi - Xj) 2 + ( Vi -yj) 2 + (Zi - Zj) 2 . 



For convenience, we take e = 1/4 and a = 1. Since Vjv only 
depends upon the pairwise distances, the set of SPs is in- 
variant under overall translation and rotation. Thus, we fix 



xi = Vi 



2/2 = Z2 = z 3 = 0. 



(1) 



Now, to construct a polynomial system equivalent to 
VVjv = 0, we a dd variables Ri t j with polynomial equations 



Ri.j ((xi - xj) 2 + (yi ~ y 3 ) 2 + (zi - Zjf) 



That is, Ri tj = l/rf^ so that V N = J2 
simplicity, we define Ri 



i<3 



(2) 



For 



SPs, we consider the polynomial system 



/jv(x, y, z,Rij) = 



Rj t i for i ^ j. Hence, for the 



E 
E 



Xi), 

l) {Vj ~ Vi), 
l) ( z j ~ z i), 



Ri,o ((xi - Xj) 2 + (yi - yj) 2 + (zi - z ) 2 ) - 1, 



i = 2, 
i = 3, 
i = 4, 
i < j 



. . , N 
..,7V 
. . , N 



An extensive search for minima and saddle points was car- 
ried out in (23} for this model up to N = 13 along with 
a corresponding search for minima and saddles of index 
one (transition states) for N = 14 in (28|. All of the 
minima and transition states, available for download at 
|http : //doye . chem . ox . ac . uk/net works/LJ ri7html| were 
obtained using numerical methods and hence they are nu- 
merical approximations. 

To certify these solutions, we first translated and rotated 
each so that condition ([T]) holds, and computed R4J based 
on ©. The downloaded points are provided to 10 deci- 
mal places, and many of them were not certifiable. We 
performed two Newton iterations using 96-bit precision to 
improve both the precision and accuracy so that ([T]) and @ 
hold. Finally, using the resulting points, we employed al- 
phaCertified to compute upper bounds on a(/jv, •)) which 
we summarize in the following table. In particular, this 
table shows that, for TV = 7, . . . , 14, each numerical ap- 
proximation of the minima and transition states does in- 
deed correspond to a certified approximate solution. For 
N = 14, the values of [log 10 7(/i4, •)] for the minima and 
transition states suggest that we should use numerical ap- 
proaches that approximate the coordinates of each station- 
ary point to at least 15 decimal places. Hence, certification 
can also provide insight into convergence conditions, which 
hitherto have been chosen based on physical intuition. 

Since we performed two Newton iterations prior to certifica- 
tion, we need to perform an a posteriori verification that we 
still have distinct solutions. This was accomplished using 
the triangle inequality, as discussed above, with the maxi- 
mum value of /3(/jvj •) and the minimum pairwise distance 
between the x, y, z coordinates of the approximations. To 
summarize, the following table proves that the numerically 
approximated SPs are indeed in the quadratic convergence 
basin of distinct SPs. 
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maximum 




maximum 


minimum 




number 


upper bound 


maximum 


upper bound 


pairwise 


N 


of points 


of a(f N , ■) 


0(fN,') 


of 7(/iV,-) 


distance 


7 


16 


6.82 ■ 10~ 19 


1.03 • 10" 28 


1.01 ■ 10 1() 


0.4354 


8 


50 


2.03 • 10" 19 


1.28 • 10" 28 


1.60 ■ 10 9 


0.4268 


9 


186 


3.55 • 10~ 17 


1.77- 10" 27 


5.46 ■ 10 1() 


0.0559 


10 


699 


2.86 • 10~ 14 


3.94 • 10~ 24 


9.87 ■ 10 10 


0.0600 


11 


2594 


6.40 ■ 10~ 15 


2.80 ■ 10" 26 


6.04 ■ 10 11 


0.0556 


12 


9122 


1.05 ■ 10" 9 


1.63 • 10" 21 


4.28 ■ 10 12 


0.0093 


13 


30265 


2.48 • 10~ 12 


2.84 • 10~ 23 


2.16 ■ 10 13 


0.0081 


14 


91415 


4.54 ■ 10~ 8 


1.50 • 10~ 19 


3.58 ■ 10 14 


0.0087 



Conclusion: Numerical approximate solutions obtained 
from standard non-linear optimization methods may lie 
in the linear convergence basin, or even in a chaotic re- 
gion, instead of the desired quadratic region of convergence. 
Hence, the numerical approximation may turn out to be a 
non-solution of the system when more Newton iterations 
are performed, which could change the scientific conclu- 
sions drastically. We have demonstrated several examples 
of such behaviour. To mitigate such problems, we shown 



how Smale's a-theory can be used to certify that a numer- 
ical approximation is in the quadratic convergence region 
of a solution, to determine if two points correspond to dis- 
tinct solutions, and to determine if the corresponding solu- 
tion is real. As a practical demonstration of the approach, 
we have refined and then certified all the known minima 
and transition states for the Lennard- Jones potential for 
up to 14 atoms. This is the first certification conducted for 
a set of physically relevant atomic structures that we are 
aware of, and it provides quantitative convergence crite- 
ria for geometry optimization. We also observe that for the 
stationary points of the Lennard- Jones potential, the size of 
the quadratic convergence basin decreases as N increases. 
All of these new insights should be applicable throughout 
molecular science and studies of soft and condensed mat- 
ter, wherever stationary points are considered to analyze 
structure, dynamics and thermodynamic properties. 
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